一文看懂线性回归【比较详细】(内含源码)

喜欢我们，点击上方AINLPer，关注一下，极品干货即刻送达！

引言

    最近作者网上看了很多关于线性回归的帖子，个人感觉比较乱，所以打算自己整理一版，希望能把线性回归说的更明白一些，为此我还整理了关于线性回归的脑图如下图所示，如果哪里有不对的地方，欢迎批评指正。

本文源码获取后台回复：线性回归源码

正文开始

1线性回归解释  1.1 线性回归定义

    线性回归（英语：linear regression）是利用称为线性回归方程的最小二乘函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。这种函数是一个或多个称为回归系数的模型参数的线性组合。只有一个自变量的情况称为简单回归，大于一个自变量情况的叫做多元回归（multivariable linear regression）。【来自维基】

1.2 线性回归的应用

    线性回归有很多实际用途。分为以下两大类：

    1、如果目标是预测或者映射，线性回归可以用来对观测数据集的和X的值拟合出一个预测模型。当完成这样一个模型以后，对于一个新增的X值，在没有给定与它相配对的y的情况下，可以用这个拟合过的模型预测出一个y值。

    2、给定一个变量y和一些变量，这些变量有可能与y相关，线性回归分析可以用来量化y与之间相关性的强度，评估出与y不相关的，并识别出哪些的子集包含了关于y的冗余信息。

    我们大多数情况下是使用第一种方式做预测，即通过已知的x，y来预测模型，进而预测未来给x，输出y的值。

1.3线性回归求解方法

    1.3.1、最小二乘法（x矩阵必须可逆）

    详情可以参考：                            https://blog.csdn.net/u014445499/article/details/107209837/

    1.3.2、梯度下降法

    对于一个函数y=f(x)，我们要求解该函数的极值点，高中的时候就学过，通过解方程f′(x)=0，就可以得到函数的极值点(x0,y0)。不过对于计算机来说，它可不会解方程。但是它可以凭借强大的计算能力，一步一步的去把函数的极值点『试』出来。

    首先要知道需要优化的目标函数为：

    其中： (x)= x= + ，对求导数得到(该公式对于多项式函数同样适用)：

    最后，梯度更新如下：

    下图就是一个最经典的梯度下降示意图。

    梯度下降法具体介绍可以参考：    https://blog.csdn.net/u013898698/article/details/120720623

    梯度下降分类

    很多视频、文章都会提到，根据数据迭代的计算方式对梯度下降分类，主要可分为以下3类：

    1）全批量梯度下降（BGD）

        1、计算量大，参数更新慢，对内存的要求很高，不能以在线的形式训练模型，也就是运行时不能加入新样本

        2、可以得到全局最优解，参数更新比较稳定，收敛方向稳定。

    2）随机梯度下降（SGD）：

        1、运算速度很快，同时能够在线学习。

        2、随机梯度下降参数更新的过程震荡很大，目标函数波动剧烈，参数更新方向有很大的波动。

        3、其较大的波动可能收敛到比批量梯度下降更小的局部极小值，因为会从一个极小值跳出来。

    3）小批量梯度下降（MBGD）：

        1、降低了更新参数的方差，使得收敛过程更加的稳定 

        2、能够利用高度优化的矩阵运算，很高效的求得每小批数据的梯度。

1.4 梯度下降和最小二乘法对比

    梯度下降法优缺点：

        1、当特征数n很大的时候，效果比较好。

    梯度下降法--缺点：

        1、需要选择学习率 alpha。

        2、需要多次的迭代。

        3、特征值取值范围比较大需要进行特征缩放（归一化）。

    最小二乘法--优点：

        1、当特征数n很大的时候，运算的比较慢，求逆矩阵的复杂度为:

    最小二乘法--缺点

        1、不需要选择学习率 alpha

        2、没有多次的迭代

        3、没有特征缩放（归一化）

1.5两种求解算法选择：

    在深度学习中，因为涉及参数特别多，动不动涉及的参数高达几十几百万，所以基本上全部都是通过梯度下降的方法求解。后面会有这种方法的代码实现，其实这块主要帮助大家加深对梯度下降的理解，在实际项目中，主要还是通过mini-batch的思想来更新梯度。

    以此衍生出来的方法也有很多，例如：SGD、SGDM（SGD+Momentum）、Adagrad、RMSProp、Adma(2015年)、swats（2017）、AMSGrad(2018)、AdaBound(2019)、Cycle-LR、OneCycle-LR，目前常用的主要有：SGDM、Adma、SWATS。这里就不展开来讨论了，后面会专门出一期方法对比的文章。

2单变量线性回归及梯度下降    本节主要介绍单变量回归模拟两种求解方法：正规方程求解、梯度下降求解，另外还针对三种梯度下降的算法进行分析对比。单变量训练数据生成

    为了模拟单变量线性回归，这里通过函数y=4.7x+2 来构造的x、y数据，作为训练数据，其中y在原有函数的基础上新增了随机变量，具体代码和生成的图形如下所示：

# 基于y=4.7x+2生成x，y训练数据
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
x=2*np.random.rand(100,1)
y=2+4.7*x+np.random.rand(100,1)
plt.plot(x,y,'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.axis([0,2,0,15])
plt.show()

最小二乘法求解

    通过回归公式来计算权重和偏置项theta，并画出求解线性函数，如红线所示。

x_b=np.c_[np.ones((100,1)),x]
theta=np.linalg.inv(x_b.T.dot(x_b)).dot(x_b.T).dot(y)
# 两个点画一条线
xx=np.array([[0],[2]])
xx_tmp=np.c_[np.ones((2,1)),xx]
yy=xx_tmp.dot(theta)
print('方程直接求解w、b的值分别为{0},{1}'.format(theta[1],theta[0]))
plt.plot(xx,yy,'r--')
plt.plot(x,y,'b.')
plt.xlabel('x')
plt.ylabel('y')
plt.axis([0,2,0,15])
plt.show()

批量梯度下降求解

    批量梯度下降主要是在整个数据集，通过对所有的样本的计算来求解梯度的方向。这里我们还设置了不同学习率（0.01、0.05、0.1、0.5）对实际训练结果的影响，相关代码如下：

theta_for_bgd=[]
# 批量梯度求解函数
def diff_lr_for_gradent(iterations,lr):
    j_save=[]
    x_b=np.c_[np.ones((100,1)),x]
    m=len(x_b)
    theta=np.random.rand(2,1)
    plt.plot(x,y,'b.')
    for iteration in range(iterations):
        y_pred=x_b.dot(theta)
        # 计算损失函数并存放到[]里面
        J=1/(2*m)*(y-y_pred).T.dot(y-y_pred)
        j_save.append(float(J[0]))   
        # 画出拟合直线
        plt.plot(x,y_pred,'r-')
        # 梯度计算及更新
        gradent=2/m*x_b.T.dot(x_b.dot(theta)-y)
        theta=theta-lr*gradent
        if lr==0.05:
            theta_for_bgd.append(theta)
    # 画出损失函数和训练次数的图像 
    plt.xlabel('X')
    plt.ylabel("Y")
    plt.axis([0,2,0,15])
    plt.title('lr={}'.format(lr)) 
    return  j_save 

iterations=100
plt.figure(figsize=(15,4))

plt.subplot(241)
j_save=diff_lr_for_gradent(iterations,0.01)
tmp=np.arange(iterations)
plt.subplot(245)
plt.plot(tmp,j_save,'g-')

plt.subplot(242)
j_save=diff_lr_for_gradent(iterations,0.05)
tmp=np.arange(iterations)
plt.subplot(246)
plt.plot(tmp,j_save,'g-')

plt.subplot(243)
j_save=diff_lr_for_gradent(iterations,0.1)
tmp=np.arange(iterations)
plt.subplot(247)
plt.plot(tmp,j_save,'g-')

plt.subplot(244)
j_save=diff_lr_for_gradent(iterations,0.5)
tmp=np.arange(iterations)
plt.subplot(248)
plt.plot(tmp,j_save,'g-')
plt.show()

    如上图，红色线是拟合线、绿色线表示了迭代次数和损失函数的关系。    可以看出来lr越小，其拟合的直线向中间逼近的比较慢（损失函数收敛慢），随着lr的增加，逼近的速度也越快（损失函数收敛快）。    但是当lr很大的时候，反而不能趋向实际要拟合的直线（损失函数会变得发散），这说明了在训练过程中要选择合适的lr。

随机梯度下降（SGD）

    随机是指每次迭代的数据是对整个数据空间的随机采样，只要采样够随机，随机梯度的期望和真实梯度是相等的。放在模型训练里，哪怕每次迭代就一个数据点，只要多迭代几个batch，最终效果一定是一样的。相关代码结果如下图所示：

theta_for_sgd=[]
x_b=np.c_[np.ones((100,1)),x]
m=len(x_b)
epoches=50 #训练次数
theta=np.random.rand(2,1)
#batch=m # 假设batch的大小和训练数据的大小一
# lr学习率更新
def lr_update(iterations):
    return 5/(50+iterations)  
for epoch in range(epoches):
    for i in range(m):
        # 在训练数据中随机的选择x，y
        data_index=np.random.randint(m)
        xi=x_b[data_index:data_index+1]
        yi=y[data_index:data_index+1]
        # 依据该条数据更新梯度
        gradent=2*xi.T.dot(xi.dot(theta)-yi)
        lr=lr_update(epoch*m+i)
        theta=theta-lr*gradent
        theta_for_sgd.append(theta)
        # 画出部分图用来展示
        if epoch<5 and i< 10 :
            y_pred=x_b.dot(theta)
            plt.plot(x,y_pred,'r-')
plt.plot(x,y,'b.')
plt.xlabel('X')
plt.ylabel('Y')
plt.axis([0,2,0,15])
plt.title('SGD')
plt.show()

    另外这里还引入了一个自适应学习率的概念，针对学习率我们希望是一开始的时候步长比较长，后来慢慢越走越小，进而达到损失函数收敛的最佳参数。如上面代码中的lr_update()函数所示。

小批量梯度下降（MBGD）

    把数据分为若干个批，按批来更新参数，利用每一批的数据来计算梯度的方向，下降起来就不容易跑偏，减少了随机性。在实际操作的时候，需要重新对训练数据做索引打乱操作（shuffle）。相关代码结果展示如下：

theta_for_mbgd=[]
epoches=50  # 50次迭代训练
batch=5 # 每次训练批次为5
m=len(x_b)
# lr学习率更新函数
def lr_update(t):
    return 5/(50+t)
t=0
# 初始化theta
theta=np.random.rand(2,1)
for epoch in range(epoches):
    # 对训练数据打乱--洗牌
    shuffle_indexes=np.random.permutation(m) 
    x_b_shuffle=x_b[shuffle_indexes]
    y_shuffle=y[shuffle_indexes]
    
    for i in range(0,m,batch):
        t=t+1     
        x_b_tmp=x_b_shuffle[i:i+batch]
        y_tmp=y_shuffle[i:i+batch]
        # 计算梯度
        gradent=2/batch*x_b_tmp.T.dot(x_b_tmp.dot(theta)-y_tmp)
        theta=theta-lr_update(t)*gradent
        theta_for_mbgd.append(theta)
        # 画出部分图像
        if epoch<5 and t< 10 :
            y_pred=x_b.dot(theta)
            plt.plot(x,y_pred,'r-')
plt.plot(x,y,'b.')
plt.axis([0,2,0,15])
plt.xlabel('X')
plt.ylabel('Y')
plt.title('MBGD')
plt.show()

三种梯度下降方法对比

    下面画出了theta在三种梯度下降方法中的变化。代码及截图如下图所示：

    BGD红色曲线，沿着曲线直接到达终点，没有出现起伏波动；

    SGD蓝色曲线，最终达到终点，但是起伏比较大；

    MBGD绿色曲线，达到终点，有起伏，但是起伏不大。

theta_for_bgd=np.array(theta_for_bgd)
theta_for_sgd=np.array(theta_for_sgd)
theta_for_mbgd=np.array(theta_for_mbgd)
plt.figure(figsize=(8,4))
plt.subplot(1,3,1)
plt.plot(theta_for_bgd[:,0],theta_for_bgd[:,1],'r-s',linewidth=1,label='bgd')
plt.legend(loc='upper left')
plt.subplot(1,3,2)
plt.plot(theta_for_sgd[:,0],theta_for_sgd[:,1],'b-*',linewidth=1,label='sgd')
plt.legend(loc='upper left')
plt.subplot(1,3,3)
plt.plot(theta_for_mbgd[:,0],theta_for_mbgd[:,1],'g-+',linewidth=1,label='mbgd')
plt.legend(loc='upper left')
plt.show()

三种梯度算法总结

全批量梯度下降（BGD）

1、计算量大，参数更新慢，对内存的要求很高，不能以在线的形式训练模型，也就是运行时不能加入新样本

2、可以得到全局最优解，参数更新比较稳定，收敛方向稳定。

随机梯度下降（SGD）：

1、运算速度很快，同时能够在线学习。

2、随机梯度下降参数更新的过程震荡很大，目标函数波动剧烈，参数更新方向有很大的波动。

3、其较大的波动可能收敛到比批量梯度下降更小的局部极小值，因为会从一个极小值跳出来。

小批量梯度下降（MBGD）：

1、降低了更新参数的方差，使得收敛过程更加的稳定

2、能够利用高度优化的矩阵运算，很高效的求得每小批数据的梯度

3多项式线性回归

    之前的时候，我们采用的直线方程：y=4.7x+2然后加上随机变量生成测试数据，这里，我们首先通过y=2x^2+x+5来生成测试数据。生成的数据图像如下图所示：

tip：实际应用要处理的数据可能还要离散、无规律，这里主要是用来解释多项式回归问题。

训练数据模拟生成

# 生成y=3x^3+2x^2+x+5数据测试数据
m=100 # x的维度
x_2=4*np.random.rand(m,1)
y_2=2*x_2**2+x_2+5+2*np.random.randn(m,1) # 后面添加了一个正态分布
plt.plot(x_2,y_2,'r.')
plt.show()

# 按照上面数据，我们还是通过一元线性函数进行求解如下:
epoches=20
batch_size=8
x_2_b=np.c_[np.ones((100,1)),x_2]
# print(x_2_b)
theta=np.random.rand(2,1)
for epoch in range(epoches):
    # 数据shuffle
    shuffle_indexes=np.random.permutation(m) 
    x_b_shuffle=x_2_b[shuffle_indexes]
    y_shuffle=y_2[shuffle_indexes]
    for i in range(0,m,batch_size):
        x_b_tmp=x_b_shuffle[i:i+batch_size]
        y_tmp=y_shuffle[i:i+batch_size]
        # 计算梯度
        gradent=2/batch_size*x_b_tmp.T.dot(x_b_tmp.dot(theta)-y_tmp)
        theta=theta-lr_update(t)*gradent
plt.plot(x_2,y_2,'r.')
y_pred=x_2_b.dot(theta)
plt.plot(x_2,y_pred,'b*')
plt.show()      

    通过上图可以发现，通过直线的方法并不能够很好的拟合出曲线，这个时候就要选择多元函数进行拟合，这就是前面介绍的多项式回归,即利用y=ax^2+bx+c进行拟合，相关代码和结果图如下所示：

epoches=500
batch_size=16
x_2_c=np.c_[np.ones((100,1)),x_2,x_2**2] # 这里可以通过添加x_2**3\x_2**4\....来计算不同degree的结果
# print(x_2_b)
theta=np.random.rand(3,1)
for epoch in range(epoches):
    # 数据shuffle
    shuffle_indexes=np.random.permutation(m) 
    x_b_shuffle=x_2_c[shuffle_indexes]
    y_shuffle=y_2[shuffle_indexes]  
    for i in range(0,m,batch_size):
        x_b_tmp=x_b_shuffle[i:i+batch_size]
        y_tmp=y_shuffle[i:i+batch_size]
        # 计算梯度
        gradent=2/batch_size*x_b_tmp.T.dot(x_b_tmp.dot(theta)-y_tmp)
        theta=theta-lr_update(t)*gradent      
plt.plot(x_2,y_2,'r.')
x2y_pred=x_2_c.dot(theta)
# print(xy_pred)
plt.plot(x_2,x2y_pred,'b.')
plt.show()      

4多变量线性回归    当一个函数不仅仅和一参数有关系的时候（例如：吴恩达教程里面房价可能和房子面积、房子建设时间、房子位置等；李宏毅教程里面的神奇宝贝cp值可能和神奇宝贝的属性值、成长系数、防御值等），就涉及到了多项式回归。

    这里以函数z=3x+6y+2创造数据集，相关代码及视图如下：

# 依据z=3x+6y+2加上随机数字来生成训练数据（x,y,z），后面以该组数据作为训练数据进行训练，让训练结果无线逼近z=3x+6y+2
import matplotlib.pyplot as plt
import numpy as np
m=50
x=np.linspace(1,20,m)
y=np.linspace(1,20,m)
x,y=np.meshgrid(x,y)
x1=x.reshape(-1,1)
y1=y.reshape(-1,1)
z1=3*x1+6*y1+2+5*np.random.randn(m*m,1)
z=z1.reshape(m,m) # 此时构造的训练数据集为：x1(m*m,1),x2(m*m,1),z1((m*m,1))
# 设置图像大小
plt.figure(dpi=100)
ax=plt.axes(projection='3d')
ax.plot_surface(x,y,z)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

    依据上面生成的数据，通过梯度下降的方法训练得到的结果为：y=3.03x+6.01y+1.48,会发现这和y=3x+6y+2是非常接近的，另外下面两张图分别也展示了实际生成的平面、以及随着迭代次数增加loss的变化曲线。相关代码及视图如下所示。

epoches=100
batch_size=256
theta_recode=[]
theta=np.random.rand(3,1)
x_y_b=np.c_[np.ones((m*m,1)),x1,y1]
t=1
y_loss=[]
# lr学习率更新函数，学习率一开始选择的比较大，后面逐渐变小。
def lr_update(t):
    return 5/(1000+t)
# 损失函数
def costFunction(X,y,theta):
    inner =np.power( X @ theta - y, 2)
    return np.sum(inner) / (2 * len(X))
for epoch in range(epoches):
    shuffle_indexes=np.random.permutation(m*m) 
    x_y_b_shuffle=x_y_b[shuffle_indexes]
    z_shuffle=z1[shuffle_indexes]
    for i in range(0,m*m,batch_size):
        x_y_b_temp=x_y_b_shuffle[i:i+batch_size]
        z_temp=z_shuffle[i:i+batch_size]
        gradent=1/batch_size*x_y_b_temp.T.dot(x_y_b_temp.dot(theta)-z_temp)
        theta=theta-lr_update(t)*gradent 
        t+=1
    loss=costFunction(x_y_b,z1,theta)
    y_loss.append(loss)
print('theta={}'.format(theta))
plt.figure(dpi=150)
ax=plt.subplot(1,2,1,projection='3d')
z_end=theta[1]*x+theta[2]*y+theta[0]
# 画图对比，之前的图和现在图对比
ax.plot_surface(x,y,z)
ax.plot_surface(x,y,z_end,color='r')
plt.subplot(1,2,2)
x_loss=np.linspace(1,epoches,epoches)
plt.plot(x_loss,y_loss,'g-')
plt.subplots_adjust(wspace =0.5, hspace =0)
plt.show()

    在实际应用中，多变量线性回归可能不知两个特征x，y，可能还有更多的值，但是都可以通过类似的方法训练。

5文章代码获取    本文源码获取后台回复：线性回归源码

长按识别下方二维码关注我们

资料整理不易，帮忙点个【赞】、【在看】吧